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Abstract This chapter describes some of the building blocks to ensure a higher level of 
confidence in the predictability and reliability (PAR) of numerical simulation of 
multiscale complex nonlinear problems. The focus is on relating PAR of nu- 
merical simulations with complex nonlinear phenomena of numerics. To isolate 
sources of numerical uncertainties, the possible discrepancy between the chosen 
partial differential equation (PDE) model and the real physics and/or experimen- 
tal data is set aside. The discussion is restricted to how well numerical schemes 
can mimic the solution behavior of the underlying PDE model for finite time 
steps and grid spacings. The situation is complicated by the fact that the avail- 
able theory for the understanding of nonlinear behavior of numerics is not at 
a stage to fully analyze the nonlinear Euler and Navier-Stokes equations. The 
discussion is based on the knowledge gained for nonlinear model problems with 
known analytical solutions to identify and explain the possible sources and reme- 
dies of numerical uncertainties in practical computations. Examples relevant to 
turbulent flow computations are included. 
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liability of Numerical Simulations, Under-Resolved Grids, Spurious Numerical 
Solutions, Nonlinear Dynamics, Spurious Bifurcations. 

1. Introduction 

The last two decades have been an era when computation is ahead of analysis 
and when very large scale practical computations are increasingly used in poorly 
understood multiscale complex nonlinear physical problems and non-traditional 
fields. Ensuring a higher level of confidence in the predictability and reliability 
(PAR) of these numerical simulations could play a major role in furthering 
the design, understanding, affordability and safety of our next generation air 
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and space transportation systems, and systems for planetary and atmospheric 
sciences, and astrobiology research. In particular, it plays a major role in the 
success of the US Accelerated Strategic Computing Initiative (ASCI) and its 
five Academic Strategic Alliance Program (ASAP) centers. Stochasticity stands 
alongside nonlinearity and the presence of multiscale physical processes as a 
predominant feature of the scope of this research. The need to guarantee PAR 
becomes acute when computations offer the ONLY way of generating this type 
of data limited simulations, the experimental means being unfeasible for any of 
a number of possible reasons. Examples of this type of data limited problem 
are: 


■ Stability behavior of re-entry vehicles at high speeds and flow conditions 
beyond the operating ranges of existing wind tunnels 

■ Flow field in thermo-chemical nonequilibrium around space vehicles 
traveling at hypersonic velocities through the atmosphere (lack sufficient 
experimental or analytic validation) 

■ Aerodynamics of aircraft in time-dependent maneuvers at high angles 
of attack (free of interference from support structures, wind-tunnel walls 
etc., and able to treat flight at extreme and unsafe operating conditions) 

■ Stability issues of unsteady separated flows in the absence of all the un- 
wanted disturbances typical of wind-tunnel experiments (e.g., geometri- 
cally imperfect free-stream turbulence) 

This chapter describes some of the building blocks to ensure a higher level of 
confidence in the PAR of numerical simulation of the aforementioned multiscale 
complex nonlinear problems, especially the related turbulence flow computa- 
tions. To isolate the source of numerical uncertainties, the possible discrepancy 
between the chosen model and the real physics and/or experimental data is set 
aside for the moment. We concentrate only on how well numerical schemes can 
mimic the solution behavior of the underlying partial different equations (PDEs) 
for finite time steps and grid spacings. Even with this restriction, the study of 
PAR encompasses elements and factors far beyond what is discussed here. It is 
important to have a very clear distinction of numerical uncertainties from each 
source. These include but are not limited to (a) stability and well-posedness 
of the governing PDEs, (b) type, order of accuracy, nonlinear stability, and 
convergence of finite discretizations, (c) limits and barriers of existing finite 
discretizations for highly nonlinear stiff problems with source terms and forc- 
ing, and/or for long time wave propagation phenomena, (d) numerical boundary 
condition (BC) treatments, (e) finite representation of infinite domains, (f) solu- 
tion strategies in solving the nonlinear discretized equations, (g) procedures for 
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obtaining the steady-state numerical solutions, (h) grid quality and grid adapta- 
tions, (i) multigrids, and (j) domain decomposition (zonal or multicomponent 
approach) in solving large problems. At present, some of the numerical uncer- 
tainties can be explained and minimized by traditional numerical analysis and 
standard CFD practices. However, such practices, usually based on linearized 
analysis, might not be sufficient for strongly nonlinear and/or stiff problems. 
We need a good understanding of the nonlinear behavior of numerical schemes 
being used as an integral part of code verification, validation and certification. 

A major stumbling block in genuinely nonlinear studies is that unlike the 
linear model equations used for conventional stability and accuracy consider- 
ations in time-dependent PDEs, there is no equivalent unique nonlinear model 
equation for nonlinear hyperbolic and parabolic PDEs for fluid dynamics. On 
one hand, a numerical method behaving in a certain way for a particular nonlin- 
ear differential equation (DE) (PDE or ordinary differential equation (ODE)) 
might exhibit a different behavior for a different nonlinear DE even though 
the DEs are of the same type. On the other hand, even for simple nonlinear 
model DEs with known solutions, the discretized counterparts can be extremely 
complex, depending on the numerical methods. Except in special cases, there 
is no general theory at the present time to characterize the various nonlinear 
behaviors of the underlying discretized counterparts. Herein, the discussion 
is based on the knowledge gained for nonlinear model problems with known 
analytical solutions to identify and explain the possible sources and remedies 
of numerical uncertainties in practical computations. 

The term “discretized counterparts” is used to mean the finite difference 
equations resulting from finite discretizations of the underlying DEs. Here “dy- 
namics” is used loosely to mean the dynamical behavior of nonlinear dynamical 
systems (continuum or discrete) and “numerics” is used loosely to mean the 
numerical methods and procedures in solving dynamical systems. We empha- 
size here that in the study of the dynamics of numerics, unless otherwise stated, 
we always assume the continuum (governing equations) is nonlinear. 

Outline: Section 2 discusses the sources of nonlinearities and the knowledge 
gained from studying the dynamics of numerics for nonlinear model problems. 
Sections 3-5 discuss some of the relevant issues and building blocks for a more 
reliable (and predictability) numerical simulation in more details. Section 6 
shows examples of spurious numerics relevant to turbulent flow computations. 
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2. Sources of Nonlinearities and Knowledge Gained from 
Nonlinear Model Problems 

Two of the building blocks for the PAR of numerical simulations are to 
identify all the sources of nonlinearities and to isolate the elements and issues 
of numerical uncertainties due to these nonlinearities. 

Sources of Nonlinearities: The sources of nonlinearities that are well known 
in computational fluid dynamics (CFD) are due to the physics. Examples of 
nonlinearities due to the physics are convection, diffusion, forcing, turbulence 
source terms, reacting flows, combustion related problems, or any combination 
of the above. The less familiar sources of nonlinearities are due to the numerics. 
There are generally three major sources: 

■ Nonlinearities due to time discretizations - the discretized counterpart 
is nonlinear in the time step. Examples of this type are Runge-Kutta 
methods. If fixed time steps are used, spurious steady-state or spurious 
asymptotic numerical solutions can occur, depending on the the initial 
condition (IC). Linear multistep methods (LMMs) (Butcher 1987) are 
linear in the time step, and they do not exhibit spurious steady states. See 
Yee & Sweby (1991-1997) and references cited therein for the dynamics 
of numerics of standard time discretizations. 

■ Nonlinearities due to spatial discretizations - in this case, the discretized 
counterpart can be nonlinear in the grid spacing and/or the scheme. Ex- 
amples of nonlinear schemes are the total variation diminishing (TVD), 
essentially nonoscillatory (ENO) and weighted ENO (WENO) schemes. 
The resulting discretized counterparts are nonlinear (in the dependent 
variables) even though the governing equation is linear. See Yee (1989) 
and references cited therein for the forms of these schemes. 

■ Nonlinearities due to complex geometries, boundary interfaces, grid gen- 
eration, grid refinements and grid adaptations (Yee & Sweby 1995) - each 
of these procedures can introduce nonlinearities even though the govern- 
ing equation is linear. 

Continuous and Discrete Dynamical Systems: Before analyzing the dynam- 
ics of numerics, it is necessary to analyze (or understand) as much as possible 
the dynamical behavior of the governing equations and/or the physical problems 
using theories of DEs, dynamical systems of DEs, and also physical guidelines. 
For stability and well-posedness considerations, whenever it is possible, it is 
also necessary to condition (not pre-condition) the governing PDEs before the 
application of the appropriate scheme (Yee & Sjogreen 200 1 a, b). The dis- 
cretized counterparts are dynamical systems on their own. They have their own 
dynamics, and they are different from one numerical method to another in space 
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and time and are different from the underlying governing PDE (Yee & Sweby 
1 993). The procedures of solving the nonlinear algebraic systems resulting from 
using implicit methods can interfere with or superpose unwanted behavior on 
the underlying scheme. Also, the same scheme can exhibit different spurious 
behavior when used for time-accurate vs. time marching to the steady states. 
For a combination of initial condition and time step, a super-stable scheme can 
stabilize unstable physical (analytic) steady states (Yee & Sweby 1993-1996). 
Super-stable scheme here refers to the region of numerical stability enclosing 
the physical instability of the true solution of the governing equation. Yee et al. 
and Yee & Sweby (1991-1997) divide their studies into two categories, steady 
state and time accurate computations. Within each category they further di- 
vide the governing PDEs into homogeneous and nonhomogeneous (i.e., with or 
without source terms), and rapidly/slowly developing and long time integration 
problems. 

Knowledge Gained from Nonlinear Model Problems. With the aid of ele- 
mentary examples, Yee et al., Yee & Sweby (1991-1997), Sweby & Yee (1994- 
1995) and Griffiths et al. (1992a, b) discuss the fundamentals of spurious be- 
havior of commonly used time and spatial discretizations in CFD. Details of 
these examples can be found in their earlier papers. These examples consist 
of nonlinear model ODEs and PDEs with known analytical solutions (the most 
straight forward way of being sure what is “really” happening with the numer- 
ics). They illustrate the danger of employing fixed (constant) time steps and 
grid spacings. They were selected to illustrate the following different nonlinear 
behavior of numerical methods: 

■ Occurrence of stable and unstable spurious asymptotes above the lin- 
earized stability limit of the scheme (for constant time steps) 

■ Occurrence of stable and unstable spurious steady states below the lin- 
earized stability limit of the scheme (for constant time steps) 

■ Stabilization of unstable steady states by implicit and semi-implicit 
methods 

■ Interplay of initial data and time steps on the occurrence of spurious 
asymptotes 

■ Interference with the dynamics of the underlying implicit scheme by 
procedures in solving the nonlinear algebraic equations (resulting from 
implicit discretizations of the continuum equations) 

■ Dynamics of the linearized implicit Euler scheme solving the time-dependent 
equations to obtain steady states vs. Newton’s method solving the steady 
equation 
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■ Spurious dynamics independently introduced by spatial and time dis- 
cretizations 

■ Convergence problems and spurious behavior of high-resolution shock- 
capturing methods 

■ Numerically induced & suppressed (spurious) chaos, and numerically 
induced chaotic transients 

■ Spurious dynamics generated by grid adaptations 

Here “spurious numerical solutions (and asymptotes)” is used to mean numer- 
ical solutions (asymptotes) that are solutions (asymptotes) of the discretized 
counterparts but are not solutions (asymptotes) of the underlying DEs. Asymp- 
totic solutions here include steady-state solutions, periodic solutions, limit cy- 
cles, chaos and strange attractors. See Thompson & Stewart (1986) and Hop- 
pensteadt (1993) for the definition of chaos and strange attractors. 

3, Minimization of Spurious Steady State via Bifurcation 
Theory 

The use of time-marching approaches to obtain steady-state numerical solu- 
tions has been considered the method of choice in computational physics for 
three decades since the pioneering work of Moretti & Abbett (1966). Moretti 
and Abbett used this approach to solve the inviscid supersonic flow over a blunt 
body without resorting to solving the steady form of PDEs of the mixed type. 
Much success was achieved in computing a variety of weakly and moderately 
nonlinear fluid flow problems. For highly complex nonlinear problems, the situ- 
ation is more complicated. The following isolates some of the key elements and 
issues of numerical uncertainties in time-marching to the steady state. Studies 
in Yee et al. (1991-1997) indicate that each of the following can affect not only 
the convergence rate but also spurious numerics other than standard stability 
and accuracy linearized numerical analysis problems. 

■ Solving an initial boundary value problem (IB VP) with unknown initial 
data 

■ Reliability of residual test 

■ Methods used to accelerate the convergence process 

■ Precondition (not condition) the governing PDE ( might introduce addi- 
tional spurious solutions beyond the solution of the underlying PDE ) 

■ Precondition the discretized counterparts ( might introduce additional 
spurious solutions beyond the underlying discretized system) 
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■ Methods in solving the nonlinear algebraic equations from implicit meth- 
ods 

■ Mismatch in implicit schemes 

■ Nonlinear schemes 

■ Schemes that are linear vs. nonlinear in time step 

■ Adaptive time step based on local error control 

It is a standard practice in time-marching to the steady-state numerical so- 
lutions to use “local time step” (varied from grid point to grid point using the 
same CFL) for nonuniform grids. However, except in finite element methods, an 
adaptive time step based on local error control is rarely used. An adaptive time 
step is built in for standard ODE solver computer packages. It enjoyed much 
success in controlling accuracy and stability for transient (time-accurate) com- 
putations. The issue is to what extent this adaptive local error control confers 
global properties in long time integration of time-dependent PDEs and whether 
one can construct a similar error control that has guaranteed and rapid con- 
vergence to the correct steady-state numerical solutions in the time-marching 
approaches for time-dependent PDEs. 

One can see that the construction of adaptive time integrators for time- 
marching to the steady states demands new concepts and guidelines and is 
distinctively different than for the time-accurate case. Straightforward applica- 
tion of adaptive time integrators for time-accurate computations might be inap- 
propriate and/or extremely inefficient for time-marching to the steady state. For 
example, an adaptive time step based on local error control for accuracy con- 
siderations is irrelevant before a steady state has been reached. Moreover, this 
type of local error control might hinder the speed of the convergence process 
with no guarantee of leading to the correct steady state. 

The twin requirements of guaranteed and rapid convergence to the correct 
steady-state numerical solution are most often conflicting, and require a full 
understanding of the global nonlinear behavior of the numerical scheme as a 
function of the discretized parameters, grid adaptation parameters, initial data 
and boundary conditions. We believe tools from bifurcation theory can help to 
minimize spurious steady-state numerical solutions. 

In many fluid problems the solution behavior is well known for certain values 
of the physical parameters but unknown for other values. For these other values 
of the parameters, the problem might become very stiff and/or strongly nonlin- 
ear, making the available numerical schemes (or the scheme in use) intractable. 
In this situation, continuation methods in bifurcation theory can become very 
useful. If possible, one should start with the physical parameter of a known or 
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reliable steady state (e.g., flow behavior is usually known for low angles of at- 
tack but not for high angles of attack). One can then use a continuation method 
such as the improved pseudo arclength continuation method of Keller (1977) 
(or the recent developments in this area) to solve for the bifurcation curve as a 
function of the physical parameter. See e.g., Doedel (2000), Shroff & Keller 
(1993) and Davidson (1997). The equations used are the discretized counter- 
part of the steady PDEs or the time-dependent PDEs. See Stephen & Shubin 
(1981) and Shubin et al. (1981) for earlier related work. If time-marching 
approaches are used, a reliable steady-state numerical solution (as a starting 
value on the correct branch of the bifurcation curve for a particular value of the 
physical parameter) is assumed. This starting steady-state numerical solution 
is assumed to have the proper time step and initial data combination and to 
have a grid spacing fine enough to resolve the flow feature. The continuation 
method will produce a continuous spectrum of the numerical solutions as the 
underlying physical parameter is varied until it arrives at a critical value p c 
such that it either experiences a bifurcation point or fails to converge. Since 
we started on the correct branch of the bifurcation curve, the solution obtained 
before that p c should be more reliable than if one starts with the physical pa- 
rameter in question with unknown initial data and tries to stretch the limitation 
of the scheme. Note that by starting a reliable solution on the correct branch of 
the bifurcation curve, the dependence of the numerical solution on the initial 
data associated with time-marching methods can be avoided before a spurious 
bifurcation occurs. 

Finally, when one is not sure of the numerical solution, the continuation 
method can be used to double check it. This approach can even reveal the true 
limitations of the existing scheme. In other words, the approach can reveal the 
critical physical parameter for which the numerical method breaks down. On 
the other hand, if one wants to find out the largest possible time step and/or 
grid spacing that one can use for a particular problem and physical parameter, 
one can also use continuation methods to trace out the bifurcation curve as a 
function of the time step and/or grid spacing. In this case, one can start with a 
small time step and/or grid spacing with the correct steady state and observe the 
critical discretized parameter as it undergoes instability or spurious bifurcation. 
Of course, this method for minimizing spurious steady states still can suffer 
from spurious behavior due to an under-resolved grid because of limited com- 
puter resources for complex practical problems. Practical guidelines to avoid 
under-resolved grids are yet another important building block toward reliable 
numerical simulations. The efficient treatment of solving the extremely large 
set of eigenvalue problems to study the type and stability of bifurcation points 
is yet another challenge. See, e.g., Fortin et al. (1996), Davidson (1997) and 
Shroff & Keller (1993) for some discussions. Consequently, further develop- 
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ment in numerical bifurcation analysis and new concepts in adaptive methods 
for time-marching to steady state hold a key to the minimization of spurious 
numerics. 

4. Source Term Treatments in Reacting Flows 

In the modeling of problems containing finite-rate chemistry or combustion, 
often, a wide range of space and time scales is present due to the reacting terms, 
over and above the different scales associated with turbulent flows, leading 
to additional numerical difficulties. This stems mainly from the fact that the 
majority of widely used numerical algorithms in reacting flows were originally 
designed to solve non-reacting fluid flow problems. Fundamental studies on 
the behavior of these schemes for reacting model problems by the author and 
collaborators were reported in Yee & Sweby (1997) and references cited therein. 
In a majority of these studies, theory from dynamical systems was used to gain 
a better understanding of the nonlinear effects on the performance of these 
schemes. The main findings are: 

■ It was shown in LeVeque and Yee (1990) that, for stiff reactions con- 
taining shock waves, it is possible to obtain stable solutions that look 
reasonable and yet are completely wrong, because the discontinuities are 
in the wrong locations. Stiff reaction waves move at nonphysical wave 
speeds, often at the rate of one grid cell per time step, regardless of their 
proper speed. There exist several methods that can overcome this diffi- 
culty for a single reaction term. For more than a single reacting term in 
fully coupled nonlinear systems, more research is needed. One impracti- 
cal way of minimizing the wrong speed of propagation of discontinuities 
is to demand orders of magnitude grid size reduction compared with what 
appears to be a reasonable grid spacing in practice. 

■ It was shown in Lafon and Yee (1991, 1992) that the numerical phe- 
nomenon of incorrect propagation speeds of discontinuities may be linked 
to the existence of some stable spurious steady-state numerical solutions. 

■ It was also shown in Lafon and Yee (1991, 1992) that various ways of 
discretizing the reaction term can affect the stability of and convergence 
to the spurious numerical steady states and/or the exact steady states. 
Pointwise evaluation of the source terms appears to be the least stable. 

■ It was shown in Yee et al. (1991) and Griffiths et al. (1992a,b) that 
spurious discrete traveling waves can exist, depending on the method of 
discretizing the source term. When physical diffusion is added, it is not 
known what type of numerical difficulties will surface. 
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From the above findings we can safely conclude that understanding the 
nonlinear behavior of numerical schemes for reacting flows and the effects 
of finite-rate chemistry on turbulence is in its infancy. However, we believe 
that knowledge gained from fundamental studies is helpful to improve some of 
the numerical difficulties that were encountered in the past. 

5. Adaptive Numerical Methods 

Another important building block for PAR is adaptive numerical methods. 
This includes adaptive temporal and spatial schemes, grid adaptation as an 
integral part of the numerical solution process, and, most of all, adaptive nu- 
merical dissipation controls. Using tools from dynamical systems, Yee et al. 
(1991-1997), Yee& Sweby (1993-1997), Griffiths et al. (1992a, b) andLafon & 
Yee (1991, 1992) showed that adaptive temporal and adaptive spatial schemes 
are important in minimizing numerically induced chaos, numerically induced 
chaotic transients and the false prediction of flow instability by direct numerical 
simulation (DNS). Their studies further indicate the need in the development 
of practical adaptive temporal schemes based on error controls to minimize 
spurious numerics due to the full discretizations. In addition, the development 
of adaptive temporal and spatial schemes based on error controls to minimize 
numerical artifacts due to the full discretizations is also needed. This is due to 
the fact that adaptive temporal or adaptive spatial schemes alone will not be 
able to provide an accurate and reliable process to minimize numerical artifacts 
for time-accurate computations. Guided by the theory of nonlinear dynamics, 
Yee et al. (1997) and Yee & Sweby (1997) presented practical examples which 
illustrated the danger of using nonadaptive temporal and spatial schemes for 
studying flow instability. 

On the subject of adaptive numerical dissipation controls, it is well known 
that reliable, accurate and efficient direct numerical simulation (DNS) of tur- 
bulence in the presence of shock waves represents a significant challenge for 
numerical methods. Standard TVD, ENO, WENO and discontinuous Galerkin 
types of shock-capturing methods for the Euler equations are now routinely used 
in high speed blast wave simulations with virtually non-oscillatory, crisp reso- 
lution of discontinuities (see reference section). For the unaveraged unsteady 
compressible Navier-Stokes equations, it was observed that these schemes are 
still too dissipative for turbulence and transition predictions. On the other hand, 
hybrid schemes, where spectral and/or higher-order compact (Pade) schemes 
are switched to higher-order ENO schemes when shock waves are detected, 
have their deficiencies. One shortcoming of this type of hybridization is that 
the numerical solution might experience a non-smooth transition at the switch 
to a different type of scheme. For 2-D and 3-D complex shock wave and shear 
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surface interactions, the switch mechanism can become non-trivial and frequent 
activation of shock-capturing schemes is possible. 

The recent work of Yee et al. (1999, 2000), Sjogreen & Yee (2000, 2001), 
and Yee & Sjogreen (200 1 a, b) indicates that appropriate adaptive numerical 
dissipation control is essential to control nonlinear instability in general, and 
for long time integration, in particular. An integrated design approach on the 
construction of adaptive numerical dissipation controls can be found in Yee & 
Sjogreen (2001a, b). 

6. Spurious Numerics Relevant to Turbulent Flow 
Computations 

This section illustrates four numerical examples that exhibit spurious nu- 
merics relevant to turbulent flow computations. The first example discusses 
spurious vortices related to under-resolved grids and/or lack of appropriate nu- 
merical dissipation/filter controls. The second example discusses spurious be- 
havior of super-stable implicit time integrators. The last two examples discuss 
spurious behavior near the onset of turbulence and/or the onset of instability 
of the steady state solution. If care is not taken, spurious bifurcation of the 
discretized counterpart and/or a numerically induced chaotic transient can be 
mistaken for the onset of physical turbulence of the governing equation. These 
examples can serve to illustrate the connection between the spurious numerical 
phenomena observed in simple nonlinear models and CFD computations. 

6.1 Spurious Vortices in Under-Resolved Incompressible 
Thin Shear Layer Flow Simulations 

Brown & Minion (1995) performed a thorough study of a second-order 
Godunov-projection method and a fourth-order central difference method for 
the 2-D incompressible Navier-Stokes equations, varying the resolution of the 
computational mesh with the rest of the physical and discretization parameters 
fixed. This is a good example of isolating the cause of spurious behavior. The 
physical problem is a doubly periodic double shear layer. The shear layers are 
perturbed slightly at the initial time, which causes the shear layer to roll up 
in time into large vortical structures. For a chosen shear layer width that is 
considered to be thin and a fixed perturbation size, they compared the solution 
for four different grid sizes (64 x 64, 128 x 128, 256 x 256, 512 x 512) 
with a reference solution using a grid size of 1024 x 1024. For the 256 x 256 
grid, a spurious vortex was formed midway between the periodically repeating 
main vortex on each shear layer. The 128 x 128 solution showed three spu- 
rious vortices along the shear layer. The spurious vortex disappeared with a 
512 x 512 mesh. They also disabled the flux limiters (a strictly upwind Fromm’s 
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method), and found the behavior to be similar. A subsequent study (Minion 
& Brown 1997) using five different formulations and six different commonly 
used schemes in CFD found similar behavior. They concluded that the spurious 
vortex is the artifact of an under-resolved grid and the behavior is caused by a 
nonlinear effect. Linking this behavior with nonlinear dynamics, we interpret 
their observation as follows. For the particular grid size and time step com- 
bination, stable spurious equilibrium points were introduced by the numerics 
into a portion of the flow field while the major portion of the flow field was 
predicted correctly. In other words, the spurious vortices are the solution of the 
discretized counterpart for that particular range of grid size and time step. The 
number of stable spurious vortices is a function of the grid size. As the grid 
spacing decreases, the spurious equilibria gradually become unstable and the 
numerical solution mimics the true solution. 

Instead of merely increasing the grid size, there are situations where under- 
resolved grids can be overcome by proper control of numerical dissipa- 
tion/filters. It was shown in Fischer & Mullen that high-order spectral ele- 
ment methods (Maday & Patera 1989), coupled with filter-based dissipation, 
can remove what is believed to be an underresolution-induced spurious vortex 
numerical solution. See Fischer & Mullen or Yee & Sjogreen (200 1 a, b) for a 
discussion. Fischer & Mullen or Yee & Sjogreen (200 1 a, b) illustrate the added 
benefit of adaptive numerical dissipation/filter controls for high order or high 
resolution shock-capturing schemes. 

6.2 Stabilizing Unstable Steady States with Implicit Time 
Integrators (Poliashenko & Yee 1995, unpublished) 

This is a joint work with Maxim Poliashenko in 1995. This unpublished 
work was presented at the 10th International Conference in Finite Element 
Methods, January 5-8, 1998, Tucson, Arizona, and also has been presented at 
various invited lectures during the last four years. We consider a 2-D lid driven 
cavity (LDC) problem. The PDE governing equations are the ideal viscous 
incompressible Navier-Stokes equations of the form 

u« + (u T V)u = -Vp + ^Au, (61) 

div u = 0, 

with boundary conditions in the domain (x, y) 


U (y = a) = 1, (6.2) 

u(y = 0) = u(a; = 0) = u(x = 1) = 0. (6.3) 

Here u is a 2-D velocity vector, p is pressure and Re is the Reynolds number, a 
dimensionless parameter of the problem that describes the relationship between 
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kinematic and viscous forces in the fluid. Here a is the cavity aspect ratio. The 
velocity vector u, pressure p and the time t are normalized with Re being 
proportional to the velocity of the lid and inversely proportional to the viscosity 
of the fluid. 

Several numerical time integrators indicated that while the steady state is 
unique and stable for small Reynolds numbers, the flow becomes time-periodic 
as Re is increased to a few thousand. Poliashenko and Aidun (1995) applied 
their direct method for computations of co-dimension one bifurcations to show 
that the steady state of the LDC problem indeed loses its stability via Hopf 
bifurcation (Thompson & Stewart 1986) as the Reynolds number increases, 
giving rise to a time-periodic solution. This bifurcation can be subcritical or 
supercritical depending on the aspect ratio a. With a given spatial resolution, 
they found that the Hopf bifurcation point is supercritical for a = 0.8 at Re = 
5220, and for a = 1.0 with Re = 7760, and subcritical for a = 1.5 at Re = 
7220, and for a = 2.0 with Re = 5120. 

For the current numerical experiment, a 47 x 47 mildly clustered finite el- 
ement mesh is used to spatially discretize the incompressible Navier-Stokes 
equations. The flow solver is the finite element code FIDAP. Nine-node quadri- 
lateral elements with biquadratic interpolation functions for velocity compo- 
nents are used. The bilinear pressure interpolation functions are projected onto 
the four Gauss points inside each element. In order to reduce the number of 
nodal unknowns, a penalty approach to remove the pressure is used. The ele- 
ments are 5 times thinner at the side walls and the bottom and 7 times thinner at 
the moving lid boundary than at the cavity center. After the spatial discretiza- 
tion, and the use of the weighted residual Galerkin method and the penalty 
formulation for the pressure, we obtain 

+ K(U)U = F. (6.4) 

at 

Here U is the global vector of system unknowns of size 2 * N where N is the 
total number of non-boundary nodes. M is a block diagonal mass matrix. The 
nonlinear matrix K represents contributions from the convective and diffusive 
terms. F is a generalized force vector which includes contribution from body 
forces. 

The dynamics of two implicit predictor-corrector time integrators are studied. 
The first is a first-order implicit Euler 


M— — — + K(U n+1 )U n+1 = F n +\ 


with the explicit Euler scheme as predictor 

U p -U n 


M 


+ K(U")U n = F". 


(6.5) 


h 


( 6 . 6 ) 
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The second is the Trapezoidal rule 


TTn+l _ TTJ> 1 1 

M + ~[K(U n+1 )U n+1 + K(U P )U P ] = -[F n+1 + F p ], (6.7) 

fbfl ^ & 

with the Adams formula as predictor 

U p = U n + %[(2 + —— — )U” - — U"- 1 ], (6.8) 

2 h n —\ 1 

where h is a constant time step, h n is a variable time step and U" is an accel- 
eration vector, approximated by 

U n = — (U n - U n_1 ) - U"- 1 . (6.9) 

hn - 1 

The local time truncation error AU n+1 is computed as follows: 


ait +1 = 


jjn+l _ ^U n+1 ) P 

3(1 + TS 1 ) 


+ 0(h\). 


( 6 . 10 ) 


Afterthenorm ||AU n+1 1| is evaluated, the next time step is computed according 
to the formula 


hn+1 /ln ^||AU n+1 || 


)t / 3 ) 


( 6 . 11 ) 


where e is a truncation error tolerance. We also set an upper limit, h max , that 
restricts growth of the time step. 

We first study the dynamics of these time integrators using a fixed time 
step (/i/c= constant). Standard Newton-Raphson and quasi-Newton iterative 
methods are used to solve the nonlinear algebraic equations for this system of 
ODEs with the predicted solution as an initial guess. With a fixed time step of 
h = 0.001, both time integrators produce a periodic solution. However, as the 
time step is increased up to h = 0.01, the first time integrator is attracted to the 
steady state solution. This phenomenon of spurious stabilization of an unstable 
steady state is very typical for implicit LMM schemes. This spurious steady 
state remains stable as h increases further. 

The dynamics of the second integrator is different. The solution remains 
time periodic up to h = 0.75. For h between 0.75 and 0.8, the solution appears 
quasiperiodic, indicating the occurrence of the secondary Hopf bifurcation. 
With h increased to 1.0, the quasiperiodic oscillations become increasingly dis- 
turbed until the solution appears very irregular for h > 1.0, which is indicative 
of numerically induced chaos. With further increases in h , more complex bi- 
furcations occur with the computed solution becoming regular again at h = 1.3 
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2-D Lid Driven Cavity (LDC) (Joint work with M. Poliashenko) 
(Different Time Integrators Exhibit Distinct Spurious Bifurcations) 


Numerical Methods : (Finite Element Code - FIDAP) 

Space : Weighted residual Galerkin method + penalty approach to remove p 
Time : Implicit Euler & trapezoidal rule 

Newton-Raphson & Quasi-Newton to solve the nonlinear algebraic Eqns. 
Explicit Euler & Adams formula as initial guesses 


BCs : v(y=a) = 1 

Rest of BCs are zero 

Grid : 47X47 

Re = 5500 

Summary : 



Variable time step trapezoidal rule 

0.15 r 



0.17 - 0.05 


= 0.15 


Implicit Euler : h=0.005 - limit cycle 

h=0.01 - steady state (stabilizing unstable SS) 


0.06 

= 0*05 


Trapezoidal : h=0.75 ■ limit cycle 

h=0.8 - quasi periodic (secondary bifurcation) 
h>0.1 - chaotic behavior 
variable time step - different spurious behavior 
(Poliashenko & Aidun, 1995) 


Figure 6.1. Two different predictor-corrector implicit methods exhibit distinct spurious bi- 
furcations. Phase portraits using a variable time step predictor-corrector implicit method 
(Trapezoidal rule as corrector and third-order Adams -Bashforth as predictor ). (Left figure: 
hmax = 0.15. Right figure: hmax = 0.05J. 
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and then returning to chaotic behavior. With fully developed chaotic behavior 
and the solution being non-smooth, more and more singular eruptions occur. 
This makes it difficult for the Newton-Raphson procedure to converge. At some 
point around h — 1.5, the Newton-Raphson method fails to converge, implying 
that there is a nonlinear instability of the chosen time integrator at this time 
step. 

Similar behavior of the numerical solution is observed if the quasi-Newton 
procedure is applied instead of the standard Newton-Raphson method in solving 
the underlying nonlinear algebraic system. However, in this case quasiperiodic 
solutions tend to become irregular with smaller time steps and the transition to 
chaos occurs earlier. 

In Poliashenko & Aidun (1995) the variable time step control version of the 
Adams-Bashforth/Trapezoidal predictor-corrector scheme with truncation error 
tolerance of 0.005 and maximum time step of 0. 15 for variable time step control 
and the quasi-Newton nonlinear solver are used. They found that after Re > 
6200, a weak modulation of the oscillation envelope occurs and the solution 
becomes quasiperiodic with the indication of a secondary Hopf bifurcation. 
As Re approaches 6500, they observed a strong resonance between the two 
independent frequencies which transformed the solution from quasiperiodic to 
strictly periodic. This phenomenon is known in dynamical systems as “phase- 
lock” on the torus. This resulted in the birth of limit cycles. This limit cycle is 
observed in a fairly wide range of Re up to 6700. At Re — 6700, the numerical 
solution exhibits more complex bifurcations and transitions to weakly turbulent, 
or chaotic motion. 

As we decrease h max < 0.12, the limit cycle shown in Fig. 6.1 is replaced 
by a stable 2-D torus which remains qualitatively unchanged as h is further 
decreased and appears to be close to the “true” solution of the ODEs. The 
latter example demonstrates that a variable time step integrator with local error 
control, although more reliable, does not guarantee no spurious numerics. 

6.3 Chaotic Transients Near the Onset of Turbulence in 

Direct Numerical Simulations of Channel Flow (Keefe 
1988, Yee & Sweby 1997) 

In addition to the inherent chaotic and chaotic transient behavior in some 
physical systems, numerics can independently introduce and suppress chaos as 
well as chaotic transients. Loosely speaking, a chaotic transient behaves like 
a chaotic solution (Grebogi et al. 1983). A chaotic transient can occur in a 
continuum or a discrete dynamical system. One of the major characteristics 
of a numerically induced chaotic transient is that if one does not integrate the 
discretized equations long enough, the numerical solution has all the character- 
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istics of a chaotic solution. The required number of integration steps might be 
far beyond those found in standard CFD simulation practice before the numer- 
ical solution can get out of the chaotic transient mode. Furthermore, standard 
numerical methods, depending on the initial data, usually experience drastic 
reductions in step size and convergence rate near a bifurcation point of the 
continuum in addition to the bifurcation points due solely to the discretized 
parameters. See Yee & Sweby (1992, 1995, 1996, 1997) for a discussion. 
Consequently, a possible numerically induced chaotic transient is especially 
worrisome in direct numerical simulations of the transition from laminar to tur- 
bulent flows. Except for special situations, it is extremely difficult to bracket 
closely the physical transition point by mere DNS of the Navier-Stokes equa- 
tions. Even away from the transition point, this type of numerical simulation is 
already very CPU intensive and the convergence rate is usually rather slow. Due 
to limited computer resources, the numerical simulation can result in chaotic 
transients indistinguishable from sustained turbulence, yielding a spurious pic- 
ture of the flow for a given Reynolds number. Consequently, it casts some doubt 
on the reliability of numerically predicted transition points and chaotic flows. 
It also influences the true connection between chaos and turbulence. See also 
Moore et al. (1990). 

Assuming a known physical bifurcation or transition point, Fig. 6.2 illus- 
trates the schematic of four possible spurious bifurcations due to constant time 
steps and constant grid spacings. This section and the next (Section 6.4) illus- 
trate the occurrence of these scenarios. Section 6.4 discusses the stability of the 
steady state (as a function of the Reynolds number) of a 2-D backward facing 
step problem using direct simulations. The present section is the computation 
by Laurence Keefe performed in the late 1980s. In 1996 we made use of the 
knowledge from continuum and discrete dynamical systems theory to interpret 
his result. We identified some of the aforementioned numerical uncertainties 
in his computations. The result is reported in Yee & Sweby 1997. 

The physical problem that Keefe considered is depicted in Fig. 6.3, where 
a flow is confined between planes at y = ±1 and is driven in the ^-direction 
by a mean pressure gradient dp/dx. The flow is characterized by a Reynolds 
number Re = U^L/v , where [Too is the mean centerline velocity, L is the 
channel half-height, and v is the kinematic viscosity. Within the channel the 
flow satisfies the incompressible Navier-Stokes equations and no-slip boundary 
conditions are applied at the walls. In the particular calculations shown here 
these equations have been manipulated into velocity-vorticity form, where one 
integrates equations for the wall-normal velocity v and normal vorticity 77 , 
and recovers the other two velocity components from the incompressibility 
condition and the definition of 77 . 
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Schematic of Possible Spurious Bifurcation 

(Assume a certain physical transition; same 1C & BC) 


1. Different temporal dlscreL (same spatial dlscreL) 

(fbnddUdt) 


Physical bit. 

(trarultton pt) 



3. Different dx (same spatial & temporal dlscreL) 

(fix'd dt) 



7 l Different dt (same spatial & temporal dlscreL) 

(finddx) 


4. Different spatial dlscreL (same temporal dlscreL) 

(tb^d til Adi) 



Figure 6 . 2 . Schematic of possible spurious bifurcation for constant time steps and grid spac- 

ings. ( 1 ) Different temporal discretizations ode i, ode2, ode 3 and ode$ ( same spatial discretiza- 
tion and the same constant dt and dx). ( 2 ) Different constant time steps dt 1, dt2, dtz and 
dt4 (same temporal and spatial discretizations, and the same constant dx). ( 3 ) Different con- 
stant grid spacings dx 1, dx 2, dx 3 and dx 4 (same spatial and temporal discretizations, and the 
same constant dt). ( 4 ) Different spatial discretizations sde \, sde 2 and sde 3 (same temporal 
discretization and the same constant dt). 
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wd 2 v - 

(6.12) 

j t n = h, + Aa 2 ^, 

(6.13) 

, dv 

f + — 0) 

dy 

(6.14) 


where 


du dw du dw 


(6.15) 


h 


V 


d_(dHi dH :i \ 
dy\dx + dz ) 


+ \9i 2 dz 2 


H 2 , 


( 6 . 16 ) 


h =^1-^1. (6.17) 

a dz dx 

Here the Hi contain the nonlinear terms in the primitive form of the Navier- 
Stokes equations and the mean pressure gradient. 

The velocity increases extremely rapidly normal to the wall, and turbulent 
channel flows are essentially homogeneous in planes parallel to the wall. The 
first requires a concentration of grid points near the wall, and the second suggests 
use of a doubly periodic domain in planes parallel to the wall. A spectral 
representation of the velocity field ( u , v , w) 


(6 - l8) 

l m n 

where the Ti(y) are Chebyshev polynomials used for the spatial discretization. 
The numerical problem then becomes dependent on a and (5 in addition to 
Re. For the time discretization, mixed explicit-implicit methods are used. 
The nonlinear terms in the equations are advanced using second-order Adams- 
Bashforth or a low storage, third-order Runge-Kutta scheme (Spalart et al. 
1991), while the viscous terms are advanced by Crank-Nicholson. 

One of the central problems in studies of wall bounded shear flows is the 
determination of when a steady laminar flow becomes unstable and transitions 
to turbulence. In dynamical systems terms, the Navier-Stokes equations always 
have a fixed point solution for low enough Reynolds numbers, but for each flow 
geometry the Reynolds number at which this fixed point bifurcates needs to be 
determined. In channel flow the fixed point solution (a parabolic velocity profile 
across the channel, u(y) = (1 — y 2 )) becomes linearly unstable at Re = 5, 772 
(Orszag 1971). However, since turbulence appears in experiments at much 
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lower Reynolds numbers, it was conjectured that this bifurcation must be sub- 
critical. Subsequent numerical solution of the nonlinear stability equations 
(Herbert 1976, Ehrenstein & Koch 1991) demonstrated this to be true, show- 
ing that limit cycle solutions with amplitude e branch back to lower Reynolds 
numbers before subsequently passing through a turning point and curving back 
toward higher Reynolds numbers. Thus for Reynolds numbers just above the 
turning point the flow equations have at least four solutions: the fixed point; 
two unstable limit cycles; and a chaotic solution (experimentally observed tur- 
bulence). Determining the location of the turning point in (a, /?, e, Re) space 
is known as the minimum-critical-Reynolds-number problem, and its solution 
is by no means complete. 

One way to investigate the turning point problem is to perform DNS of 
channel flow for conditions believed to be near this critical condition. Beginning 
with a known turbulent initial condition from higher Reynolds number, one 
integrates in time at the target Reynolds number to determine whether the flow 
decays back to the fixed point or sustains itself as turbulence. Although this may 
not be the most efficient way to bracket the turning point, it has the advantage 
that the peculiar dynamics of the flow near the turning point, whether in decay 
or sustained turbulence, are observable. This yields information about the path 
along which flows become turbulent at these low Reynolds numbers. 

Unfortunately the flow dynamics are very peculiar near the turning point, and 
extremely long chaotic transients are observed in the computations that make a 
fine determination of that point all but impossible by this method. This can be 
seen in Fig. 6.3, where a time history of the turbulent energy in a channel flow 
(energy above that in the laminar flow) is plotted for a Reynolds number of 2, 1 9 1 . 
To understand the time scale of the phenomenon some experimental facts need 
to be recalled. In typical experimental investigations of channel flow, the infinite 
transverse and streamwise extent of the ideal flow are approximated by studying 
flow in high aspect ratio (10-40) rectangular ducts that typically are 50-100 duct 
heights long. If times are non-dimensionalized by the centerline mean velocity 
Uoq and the duct half height L, then statistics on turbulence are gathered by 
averaging hot-wire data over intervals AtU^/L ~ 200. In the simulations and 
figure the time scale is based on the friction velocity u T and L, where typically 
15-20 u T ~ Uoq. Thus averaging over intervals A tu r /L ~ 10 should and does 
yield stable flow statistics that compare well with experiments. The near-wall 
velocity profile, cross- channel turbulence intensities, and Reynolds and shear 
stress distribution for the A tu r /L ~ 10 interval near the end of the transient, 
delineated by the arrows in Fig. 6.3, indicate the good comparison. In each 
case they correspond well to available experimental data. Yet look at the time 
scale of the transient; it spans A tu r /L ~ 300, thirty times longer than the time 
needed to obtain stable statistics that would convince most experimentalists 
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Chaotic Transient Near the Onset of Turbulence 

(Keefe 1996) 

Re = 2191, 32 X 33 X 32 grid, mixed explicit & implicit spectral method 
163,840 steps, transient calculation length of 409.6, 40 hrs on Cray XMP 


ResU^L/v 

U M - mean centerline velocity 
L - channel half height 
v - kinematic viscosity 



Performed the "dimension & Lyapunov exponent" study 
at midpoint of later computation (Re = 3200, 429,680 steps) 


Figure 6.3. 3-D channel flow computation by Keefe, 1988. 
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that they are viewing a fully developed turbulent channel flow. This is further 
complicated by the wide variation of the transient length, depending upon both 
the grid resolution (number of modes in the spectral representation) and the 
linearly stable time step of the integration. In fact, for fixed (a, (3, Re) it is 
possible to obtain sustained turbulence for one time step, but see it rapidly 
decay to laminar flow for another, lower value of the step. 

Extended chaotic transients near bifurcation points are not an unknown phe- 
nomenon; the "meta-chaos" of the Lorenz system is but one of many known 
examples. However, the practicalities of numerical computation in fluid dynam- 
ics usually interfere with one’s ability to discern whether transient, or sustained 
turbulence, is being calculated. The computations required to obtain the tran- 
sient plot in Fig. 6.3 needed 40 hours of single processor time on a Cray XMP, 
some ten years ago. Such a small amount of expended time was only possible 
because the spatial resolution of the calculation was relatively coarse (32 x 33 
x 32), in keeping with the large scales of the phenomena expected at these flow 
conditions. Higher resolution calculations (192 x 129 x 160) (Kim et al. 1987) 
at greater Reynolds numbers typically have taken hundreds of hours (~ 250) 
to barely obtain the A tu T /L = 10 averaging interval that is so inadequate for 
detecting transients. Because such calculations are so time consuming, one typ- 
ically chooses an integration time step that is a substantial fraction of the linear 
stability limit of the algorithm, so as to maximize the calculated "flow time" for 
expended CPU time. However, it is clear from these transient results that this 
practice has some dangers when close to critical points of the underlying con- 
tinuous dynamical system. Thus it appears that just as pseudo-time integration 
to obtain steady solutions can result in spurious results, genuine time integra- 
tion can result in chaotic transients indistinguishable from sustained turbulence, 
also yielding a spurious picture of the flow for a given Reynolds number. 

6.4 Temporal & Spatial Refinement Studies of 2-D 

Incompressible Flow over a Backward-Facing Step 

The 2-D incompressible flow over a backward-facing step has been addressed 
by many authors using a wide variety of numerical methods. Figure 6.4 shows 
the flow geometry. Fluid with constant density p and viscosity p enters the up- 
stream channel of height h with a prescribed velocity profile (usually parabolic). 
After traveling a distance l, the fluid passes over a backward-facing step of height 
s and enters the downstream channel of height H = h + s. After traveling a 
distance L downstream of the step, the fluid exits the region of interest. For 
Reynolds numbers considered here, the flow separates at the comer and forms 
a recirculating region behind the step. Additional recirculating regions form on 


Building Blocks for Reliable Simulations 


23 


Backward Facing Step 

(2-D Incompressible Flow Simulations) 



Early 90’s Controversy : Transition point Reynolds # 

• Reports of sustained unsteady flow for Reynolds # 
in the range of (250, 2500) 

• Formulations 

vortex method, unsteady Eqns. in stream function form, 
unsteady Eqns. & the associated linear-stability problem, 
unsteady Eqns. in primitive variable form 

• Numerical Methods 

All of the existing schemes in the literature 

Gresho et al. (1993) : Provided an answer to the above controversy 

(the steady solution at Re=800 is stable) 

• Kaiktsis et al. (1991) - transition to turbulent flow has occurred at Re=800 

• Torczynski (1993) - the result of Kaiktsis et al. (1991) is an artifact of 

inadequate spatial resolution 

• Torczynski’s conclusion was confirmed by a subsequent study of 

Kaiktsis etal. (1996) & Fortin et al. (1996) 


Figure 6.4. Schematic of the backward facing step problem and background. 
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the upper and subsequently the lower walls of the downstream channel as the 
Reynolds number is increased. 

Results of sustained unsteady flow from various numerical simulations have 
been reported for Reynolds numbers (Re) ranging from 250 up to 2500. The 
formulations included the vortex method, unsteady equations in stream func- 
tion form, steady equations and the associated linear-stability problem, and the 
unsteady equations in primitive variable form. The numerical methods used 
cover almost all of the existing schemes in the literature. The majority of the 
numerical results are summarized in Gresho et al. (1993). The work of Gresho 
et al. was an answer to a controversy concerning the stability of the stationary 
solution at Re — 800. It was concluded by Kaiktsis et al. (1991) that transi- 
tion to turbulent flow has occurred at Re = 800. Kaiktsis et al. examined the 
long-time temporal behavior of the flow and found that the flow is steady at 
Re = 500, time-periodic at Re — 700, and chaotic at Re = 800. Gresho et al. 
did a detailed grid refinement study using four different numerical methods and 
concluded that the backward-facing step at Re = 800 is a stable steady flow. 

In addition to the study of Gresho et al., an extensive grid refinement study of 
this flow using a spectral element method was conducted in Torczynski (1993). 
The simulated geometry and the numerical method corresponds to that of Kaik- 
tsis et al. (1991). Flow was examined at Reynolds numbers of 500 and 800. 
His systematic grid refinement study was performed by varying both the ele- 
ment size and the order of the polynomial representation within the elements. 
For both Reynolds number values with the transient computations stopped at 
t = 800, it was observed that low-resolution grid cases exhibit chaotic-like tem- 
poral behavior whereas high-resolution grid cases evolve toward asymptotically 
steady flow by a monotonic decay of the transient. The resolution required to 
obtain asymptotically steady behavior is seen to increase with Reynolds num- 
ber. These results suggest that the reported transition to sustained chaotic flow 
(Kaiktsis et al., 1991) at Reynolds numbers around 700 is an artifact of inad- 
equate spatial resolution. Torczynski’s conclusion was further confirmed by 
a subsequent study of Kaiktsis et al. (1996) and Fortin et al. (1996). Fortin 
et al. employed tools from dynamical systems theory to search for the Hopf 
bifurcation point (transition point). They showed that the flow remains steady 
at least up to Re = 1600. 

Grid Refinement Study of Torczynski (1993): In Torczynski (1993), the 
Re = pu2h/n is based on upstream conditions. The variable u is the spatial 
average of the horizontal velocity u over h. The geometry is specified to match 
that of Kaiktsis et al. (1991). The upstream channel height h and step height s 
have values of h = 1 and s — 0.94231, yielding a downstream channel height 
of H = 1.94231. The comer of the step is at ( x,y ) = (1,0). The channel 
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extends a distance L = 1 upstream from the step and a distance L = 34 
downstream from the step to preclude undue influence of the finite channel 
length on the flow at Re — 800. The following conditions are applied on 
the boundaries of the computational domain: u = v = 0 on the upper and 
lower channel walls, — p + ydu/dn = 0 and dv/dn = 0 on the outflow 
boundary, and u = [tanh(t/lQ)]iiB{y) + [1 — tanh(t/16)]up(y) and v = 0 on 
the inflow boundary and the step surface. Here, ub(v) = max[0, 3y(l — y )] 
is the correct boundary condition for flow over a backward-facing step and 
up(y) = 3(1 — y)(s + y)/( 1 + s) 3 is the Poiseuille flow observed infinitely 
far downstream whenever steady flow is asymptotically obtained. The initial 
velocity field is set equal to u = up(y) and v = 0 throughout the domain. Here 
v is the vertical velocity and p is the pressure. Thus, the above combination of 
boundary and initial conditions initially allows flow through the step surface so 
that the simulations can be initialized using an exact divergence-free solution of 
the Navier-Stokes equations. Furthermore, since the inflow boundary condition 
is varied smoothly in time from Poiseuille flow to flow over a backward-facing 
step, the flow experiences an order-unity transient that is probably strong enough 
to excite sustained unsteady behavior, if that is the appropriate asymptotic state 
for the numerical solution. 

The simulations were performed using the commercial code NEKTON v2.8, 
which employs a time-accurate spectral-element method with the Uzawa for- 
mulation (NEKTON, 1991). Let D be the dimensionality. Each element has 
N d velocity nodes located at Gauss-Lobatto Legendre collocation points, some 
of which are on the element boundaries, and (N — 2) D pressure nodes located 
at Gauss Legendre collocation points, all of which are internal. Within each 
element, the velocity components and the pressure are represented by sums 
of D-dimensional products of Lagrangian-interpolant polynomials based on 
nodal values. This representation results in continuous velocity components 
but discontinuous pressure at element boundaries. Henceforth, the quantity N 
is referred to as the element order, even though the order of the polynomials 
used to represent the velocity is N — 1. NEKTON employs mixed explicit and 
implicit temporal discretizations. To avoid solving a nonlinear nonsymmetric 
system of equations at each time step, the convective term is advanced explic- 
itly in time using a third-order Adams-Bashforth scheme. All other terms are 
treated implicitly (implicit Euler for the pressure and for the viscous terms). 

Three spectral-element grids of differing resolution, denoted L (low), M 
(medium), and H (high), are employed. Figure 6.6 shows the computational 
domain and the grid distribution of the three spectral element grids in which the 
distribution of nodes within each spectral element is not shown. The L grid 
with N = 9 is identical to the grid of Kaiktsis et al. (1991). Four general classes 
of behavior are observed for the numerical solutions. First, “steady monotonic” 
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denotes evolution of the numerical solution toward an asymptotically steady 
state. Second, “steady oscillatory” denotes evolution toward an asymptotically 
steady state with a decaying oscillation superimposed on the monotonic decay. 
Third, “unsteady chaotic” denotes irregular transient behavior of the numerical 
solution that shows no indication of evolving toward steady behavior. Fourth, 
“diverge” denotes a numerical solution terminated by a floating-point exception. 
In Fig. 6.6, the first character denotes the grid resolution L , M or H , the 
first digit indicates the Reynolds number 500 or 800 and the last two digits 
indicate the order of the spectral element being used. For example, L807 
means Re = 800 using the L grid with N — 7. 

The extensive grid refinement study of Torczynski resulted in grid-independent 
steady-state numerical solutions for both Re — 500 and Re — 800. As the 
grid resolution is reduced below the level required to obtain grid independent 
solutions, chaotic-like temporal behavior occurred. The degree of grid reso- 
lution required to obtain a grid-independent solution was observed to increase 
as the Reynolds number is increased. Figure 6.5 shows the streamlines for 
for H 809 (steady solution) and L811 (spurious time-periodic solution) and the 
corresponding grids with the distribution of the nodes of the spectral elements 
shown. 

Temporal Refinement Studies Using Knowledge from Dynamical Systems 
Theory: All of Torczynski’s numerical solutions integrate to t — 800. With the 
knowledge of possible nonlinear behavior of numerical schemes such as long 
time transients before a steady state is reached, numerically induced chaotic 
transients, numerically induced or suppressed chaos, existence of spurious 
steady states and asymptotes, and the intimate relationship among initial data, 
time step and grid spacing observed in discrete dynamical systems theory, Yee 
et al. (1997) examined the Torczynski cases in more detail. 

In the Yee et al. (1997) study, in addition to grid refinement, temporal refine- 
ments are made on all of the under-resolved grid cases to determine if these cases 
sustain the same temporal behavior at a much later time or evolve into a different 
type of spurious behavior. At t = 800, cases L506, L507, L508, L509, L811, 
M807 and M808 either exhibit “unsteady chaotic” or “steady oscillatory” be- 
havior. We integrate these cases to t — 2000 to determine if a change in 
solution behavior occurs. From the phenomena observed in Keefe’s 3-D chan- 
nel flow computation and others, t = 2000 might not be long enough for a 
long time transient or long chaotic transient to die out. There is also the po- 
tential of evolving into a different type of spurious or divergent behavior at a 
much later time. However, for this study it appears that t = 2000 is suffi- 
cient. For Re = 500, we also recomputed some of these cases with a sequence 
of At that bracketed the benchmark study of Torczynski. The At values are 
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Backward Facing Step Simulations ^ 

JBy a spectral element computer code NEKTON)^ 


3 Grids (L, M & H) 


M 

H 


Streamlines & Grids 
Re=800, t=2000 (vertical expanded 4 times) 


Correct 


Re = 800, t = 2000: N = 9, dt = 0.10 
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Figure 6.5. Three different grids, and streamlines ofH 809 and L811 and their corresponding 
grids by the spectral element method. 
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0.02,0.05,0.10,0.125,0.2,0.3,0.4, and 0.5 for Re = 500. The CFL number 
for all of these cases is above 1 for At > 0.10. The reason for the investigation 
of At = 0.3, 0.4 and 0.5 is to find out, after the transients have died out, if the 
solution converges to the correct steady state for At that are a few times larger 
than 0.10. 

Fori?e = 800, we integrate L8 11 and M808 with At = 0.10 and M807 with 
At — 0.02,0.05 and 0.10 to t = 2000. Aside from integrating to t = 2000, 
five different initial data were examined for cases M807, M809 and M 811 for 
At = 0.10 to determine the influence of the initial data and the grid resolution 
on the final numerical solution. The five initial data are: 

(a) Uniform: u, v = 0 

(b) Shear layer: u = up{y) = max[ 0, 3y(l — y)]. v = 0 

(c) Solution from solving the steady Stokes equation (with no convection 
terms) 

(d) Torczynski (1993): u = up(y) = 3(1 — y)(s + y)f{ 1 + s) 3 ,v = 0 

(e) Channel flow both upstream & downstream of step: Same as (d) except 
the boundary conditions 

The boundary conditions for (a), (b), (c) and (e) were parabolic inflow and no- 
slip at walls, whereas the boundary conditions for (d) were those of Torczynski 
(1993): u = [tanh(t/ld)]uB{y) + [1 — tanh(t/16)]up(y) and v = 0. The 
CPU time required to run the above cases ranged from less than a day to several 
days on a Sparc Center 2000 using one processor. 

The chaotic-like behavior evolves into a time-periodic solution beyond t = 
800 for L506 and L507, whereas the chaotic-like behavior evolves into a time- 
periodic solution beyond t = 800 for L811 and a divergent solution for M807. 
The “steady oscillatory” case L508 slowly evolves to the correct steady state 
with an amplitude of oscillation of 10 -5 . The oscillation is not detectable 
within the plotting accuracy. The “steady oscillatory” time evolution of M808 
is similar to that of L508. The numerical solutions with “steady oscillatory” and 
“steady monotonic” behavior at early stages of the time integration are almost 
identical at later stages of the time integration. They all converge to the correct 
steady state. The initial data study at Re = 800 with At = 0.10 is summarized 
in Table 5.5 of Yee et al. (1997). It illustrates the intimate relationship between 
initial data and grid resolution. 

Figure 6.6 shows the vertical velocity time histories at (x,y) = (30,0) 
advanced to a time of t — 2000 for M807 with At = 0.02,0.05 and 0.10, 
and L811 for At — 0.10. Case M807 diverges at t = 1909.2 for At = 0.02, 
at t = 972.4 for At = 0.05, and at t — 827.77 for At — 0.10. The time 
histories for these three time steps appear to show chaotic-like behavior if one 
stops the computations at t — 800. The bottom plot of Fig. 6.6 shows the 
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Vertical Velocity Time Histories at ( 30 , 0 ) 

Re= 800 , M & L Grids 






Figure 6.6. Vertical velocity time histories for the MSOTwith time steps 0.02, 0.05, 0.10, and 
L811 with time step 0.01 for t = 2000. 
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vertical velocity time histories advanced to a time of t — 2000 for L811 with 
At = 0.10. It shows the definite time-periodic spurious solution pattern. On 
the other hand, the time history for this case appears to show an aperiodic-like 
pattern if one stops the computation at t — 800. Note that the L809 grid case 
was used by Kaiktsis et al. (1991) and they concluded that “2-D transition” has 
already occurred at Re = 800. 

In summary, without the temporal refinement study (longer time integration), 
the L506, L507, L811 and M807 cases can be mistaken to be chaotic-like (or 
aperiodic-like) flow. Although the time history up to t = 800 appears chaotic- 
like, one cannot conclude it is chaotic without longer transient computations. 
One can conclude that with transient computations that are 2.5 times longer 
than Torczynski’s original computations, what appeared to be aperiodic-like or 
chaotic-like behavior at earlier times evolved toward either a time-periodic or 
divergent solution at later times. These temporal behaviors appear to be long 
time aperiodic-like transients or numerically induced chaotic-like transients. 
For Re — 800, five different initial data were examined to determine if the 
flow exhibits strong dependence on initial data and grid resolution. Results 
showed that the numerical solutions are sensitive to these five initial data. Note 
that the results presented pertain to the characteristic of the studied scheme 
and the direct simulations. However, if one is certain that Re — 800 is a stable 
steady flow, a non-time-accurate method such as time-marching to obtaining the 
steady-state numerical solution would be a more efficient numerical procedure. 

Spurious Bifurcation by Different Time Integrators (Henderson & Yee 1998, 
unpublished): This is a joint work with Ronald Henderson in 1998. The un- 
published work was presented at the 10th International Conference in Finite 
Element Methods, January 5-8, 1998, Tucson, Arizona, and also has been pre- 
sented at various invited lectures during the last four years. Our joint work 
illustrates the situation where solving the nonlinear terms of the Navier-Stokes 
equations by two different explicit time integrators (same implicit time integra- 
tor for the linear terms) results in spurious bifurcation. This spurious bifurcation 
is shown in Fig. 6.7 as a function of the Reynolds number. These computations 
use the implicit Euler time integrator for the linear terms. Also the same spatial 
discretization L809 is used with a fixed time step of t — 0.10. The two explicit 
time integrators are the third-order Adams-Bashforth (AB3) and a second-order 
explicit stiffly stable method (SS2) (Henderson 1999). The AB3 method ex- 
periences a spurious bifurcation near Re = 720, whereas the SS2 method 
experiences a spurious bifurcation at a larger Reynolds number near Re = 800. 
The method and the scaling for this figure can be found in Henderson (1999). 
Finding the exact location of these spurious bifurcation points requires more 
complicated computation which is not performed here. In addition, the exact 
representation of this bifurcation plot is rather complicated to explain, and is not 
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Backward Facing Step (Joint work with Henderson & Torczynski) 
(Different Time Integrators Exhibit Distinct Spurious Bifurcations) 

Numerical Methods : 

Space : Spectral element (L809 grid) 

Time : Implicit Euler (linear term) 

AD3 & SS2 (nonlinear term) 


BC_& 1C : Same as Torczynski + continuation method Spurious Bifurcations using AD3 & SS2 

(Fixed dt,L809 grid) 


Summary : (Explicit time integrators) 

* Different time integrators exhibit distinct 
spurious bifurcations (same grid, same dt) 

• As the grid size increases, spurious 
bifurcations occur at larger Re (same dt) 

Adaptive Grid Refinement : (Explicit time integrators) 
Remains laminar at least up to Re=1800 


Future Work : (Uncond. stable implicit time integrators) 

• Spurious bifurcations occur at larger Re 
than explicit methods 

• Uncond. stable implicit LMMs can stabilize 
unstable SS as dt increases (same grid) 

• Adaptive time & spatial discretizations 



Figure 6.7. Different explicit time integrators (for the nonlinear terms of the Navier-Stokes 
equations) exhibit distinct spurious bifurcation using the same spatial discretization L809 and 
a time step of 0.10. 
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important for the current discussion. They are not the main illustration for this 
study. When an adaptive version of the spectral element method (Henderson, 
1999) is used, the problem remains laminar up to Re = 1800. Future work 
which is indicated on Fig. 6.7 is planned. 

Minimization of Spurious Bifurcation by a Suitable Filter (Fischer 200 1 , un- 
published): Recently, Fischer (2001) computed the same L811 spectral element 
grid using a time integrator based on the operator integration-factor splitting 
(OlFS) developed by Maday, Patera and Rpnquist (1990). This scheme decou- 
ples the convective step from the Stokes update, thereby allowing CFL numbers 
in excess of unity. At the end of each step, Fischer applies a filter to the velocity 
that effectively scales the Nth-order Legendre modes within each element by 
(1 — a), where, typically, 0.05 < a < 0.30 (Fischer & Mullen 2001). Because 
the filter is applied on each step, its strength is a function of At as well as a. 
The spurious behavior observed by Kaiktsis et al. (1991) is cured by the filter 
and a stable steady-state numerical solution is obtained without further grid 
refinement. Figure 6.8 illustrates the velocity time histories at (30,0) by the 
filtered and un-filtered spectral element methods with At — 0.10. 

In summary. Sections 6.3 and 6.4 illustrate all of the possible scenarios of 
spurious bifurcations indicated on the schematic diagram of Fig. 6.2. The last 
scenario, discussed briefly at the beginning of Section 6.4, is quite common and 
is not shown here. See Gresho et al. (1993) and references cited therein for 
some examples. 


VI. Concluding Remarks 

Some building blocks to ensure a higher level of confidence in PAR of nu- 
merical simulations have been discussed. The discussion concentrates only on 
how well numerical schemes can mimic the solution behavior of the underlying 
PDEs. The possible discrepancy between the chosen model and the real physics 
and/or experimental data is set aside. These building blocks are based largely 
on the author’s view, background and integrated experience in computational 
physics, numerical analysis and the dynamics of numerics. They also represent 
the end result of the various studies with the author’s collaborators indicated 
in the acknowledgment Section. Among many other important building blocks 
for the PAR of numerical simulations, the author believes the following five 
building blocks are essential. The first building block is to analyze as much as 
possible the dynamical behavior of the governing equation. For stability and 
well-posedness considerations, whenever it is possible, it is also necessary to 
condition (not pre-condition) the governing PDEs before the application of the 
appropriate scheme (Yee & Sjogreen 2001a, b). The second building block is to 
understand the nonlinear behavior, limits and barriers, and to isolate spurious 
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Figure 6.8. Comparison of the filtered with the un-filtered solutions of L811 with time step 
dt = 0.10. Top: no filter; Middle: filter (filter strength 0.05); Bottom: filter (filter strength 0.3). 
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behavior of existing numerical schemes. A third building block is to include 
nonlinear dynamics and bifurcation theories as an integral part of the numerical 
process whenever it is possible. A fourth building block is to construct ap- 
propriate adaptive spatial and temporal discretizations that are suitable for the 
underlying governing equation. The last building block is to construct appro- 
priate adaptive numerical dissipation/filter controls for long time integrations, 
complex high speed turbulent and combustion simulations (Sjogreen & Yee 
2001, Yee & Sjogreen 2001). 

The need for the study of dynamics of numerics is prompted by the fact 
that the type of problem studied using CFD has changed dramatically over 
the past two decades. CFD is also undergoing an important transition, and it 
is increasingly used in nontraditional areas. But even within its field, many 
algorithms widely used in practical CFD applications were originally designed 
for much simpler problems, such as perfect or ideal gas flows. As can be 
seen in the literature, a straightforward application of these numerical methods 
to high speed flows, nonequilibrium flows, advanced turbulence modeling or 
combustion related problems can lead to wrong results, slow convergence, or 
even nonconvergent solutions. The need for new algorithms and/or modification 
and improvement to existing numerical methods in order to deal with emerging 
disciplines is evident. We believe the nonlinear dynamical approach for CFD 
can contribute to the success of this goal. 

We have revealed some of the causes of spurious phenomena due to the nu- 
merics in an attempt to improve the understanding of the effects of numerical 
uncertainties in CFD. We have shown that guidelines developed using lineariza- 
tion methods are not always valid for nonlinear problems. We have gained an 
improved understanding of long time behavior of nonlinear problems and non- 
linear stability, convergence, and reliability of time-marching approaches. We 
have learned that numerics can introduce and suppress chaos and can also in- 
troduce chaotic transients. The danger of relying on DNS to bracket closely 
the the onset of turbulence and chaos is evident. 


We illustrated with practical CFD examples that exhibit properties and quali- 
tative behavior similar to those of elementary examples in which the full dynam- 
ical behavior of the numerics can be analyzed. The observed spurious behavior 
related to under-resolved grid cases is particularly relevant to DNS and large 
Eddy simulation (LES). Spatial resolutions in DNS and LES are largely dictated 
by the computer power. These studies serve to point out the various possible 
dangers of misinterpreting numerical simulations of realistic complex flows that 
are constrained by available computing power. 


As can be seen, recent advances in dynamics of numerics show the advan- 
tage of adaptive step size error control for long time integration of nonlinear 


Building Blocks for Reliable Simulations 


35 


ODEs. Although much research is needed to construct suitable yet practical 
similar adaptive methods for PDEs, these early developments lead our way to 
future theories. There remains the challenge of constructing practical spatial 
and temporal adaptive methods for time-accurate computations, and construct- 
ing adaptive step size control methods that are suitable yet practical for time 
marching to the steady state for aerospace CFD applications. Another even 
more challenging area is the quest for an adaptive numerical scheme that leads 
to guaranteed and rapid convergence to the correct steady-state numerical so- 
lutions. 

Acknowledgments 

The author wishes to thank her collaborators Peter Sweby, Laurence Keefe, 
John Torczynski, Ronald Henderson, Maxim Poliashenko, David Griffiths, An- 
dre Lafon and Andrew Stuart, for contributing to her earlier work. The contri- 
butions of Section 6.2 by Laurence Keefe, Section 6.3 by Maxim Poliashenko, 
and Section 6.4 by Paul Fischer, John Torczynski and Ronald Henderson are 
gratefully acknowledged. Special thanks to Dennis Jespersen, Unmeel Mehta, 
Murray Tobak and Marcel Vinokur for their suggestions and critical review of 
the manuscript. 

References 

Brown, D.L. and Minion, M.L. (1995), “Performance of Under-resolved 
Two-Dimensional Incompressible Flow Simulations,” J. Comput. Phys., Vol. 
122, pp. 165-183. 

Butcher, J.C. (1987), Numerical Analysis of Ordinary Differential Equations , 
John Wiley & Son, Chichester. 

Davidson, B. (1997), Large Scale Continuation and Numerical Bifurcation 
for PDE’s, SIAM J. of Numer. Analy., Vol. 34, No. 5, pp. 2008-2027. 

Doedel, E. (2000), AUTO: Software for Continuation and Bifurcation Prob- 
lems in Ordinary Differential Equations, Concordia University, Montreal, 
Canada and Cal. Tech., Pasadena, Calif. 

Ehrenstein, U. and Koch, W. (1991), "Nonlinear bifurcation study of plane 
Poiseuille flow," J. Fluid Mech., Vol. 228, pp. 111-148. 

Fischer, P.F. and Mullen, J.S., (2001) “Filter-Based Stabilization of Spectral 
Element Methods”, Comptes Rendus de V Academie des sciences Paris, t. 332, 
Serie I - Analyse numerique, 265-270 (2001). 

Fischer, P.F., (2001) Private communication. 

Fortin, A., Jardak, M., Gervais, J.J. and Pierre, R. (1996), “Localization of 
Hopf Bifurcations in Fluid Flow Problems,” Intern. J. Numer. Meth. Fluids. 

Gresho, P.M., Gartling, D.K., Torczynski, J.R., Cliffe, K.A., Winters, K.H., 
Garrett, T.J., Spencer, A. and Goodrich, J.W. (1993), “Is the Steady Viscous In- 


36 


compressible Two-Dimensional Flow Over a Backward-Facing Step at Re=800 
Stable?” Intern. J. Numer. Meth. Fluids, Vol. 17, pp. 501-541. 

Grebogi, C., Ott, E. and Yorke, J.A. (1983), “Crises, Sudden Changes in 
Chaotic Attractors, and Transient Chaos,” Physica 7D, pp. 181-200. 

Griffiths, D.F., Sweby, P.K. and Yee, H.C. (1992a), “On Spurious Asymp- 
totes Numerical Solutions of Explicit Runge-Kutta Schemes,” IMA J. Numer. 
Anal., Vol. 12, pp. 319-338. 

Griffiths, D.F., Stuart, A.M. and Yee, H.C. (1992b), “Numerical Wave Prop- 
agation in Hyperbolic Problems with Nonlinear Source Terms,” SIAM J. of 
Numer. Analy., Vol. 29, pp. 1244-1260. 

Henderson, R. D., (1999) “Adaptive Spectral Element Methods for Turbu- 
lence and Transition, in High-Order Methods for Computational Physics”, T. J. 
Barth and H. Deconinck, editors. Springer. 

Herbert, Th. (1976), Lecture Notes in Physics 59, Springer-Verlag, p235. 

Hoppensteadt, F.C. (1993), Analysis and Simulation of Chaotic Systems, 
Springer-Verlag, New York. 

Kaiktsis, L., Kamiadakis, G.E. and Orszag, S.A. (1991), “Onset of Three- 
Dimensionality, Equilibria, and Early Transition in Flow Over a Backward- 
Facing Step,” J. Fluid Mech., Vol. 231, pp. 501-528. 

Kaiktsis, L., Kamiadakis, G.E. and Orszag, S.A. (1996), “Unsteadiness and 
Convective Instabilities in Two-Dimensional Flow Over a Backward-Facing 
Step,” J. Fluid Mech., Vol. 321, pp. 157-187. 

Keefe, L., Moin, P. and Kim, J. (1992), “The Dimension of Attractors Un- 
derlying Periodic Turbulent Poiseuille Flow,” J. Fluid Mech., Vol. 242, pp. 
1-29. 

Keefe, L. (1988-1996), unpublished; private communication. 

Kim, J., Moin, P. and Moser, R. (1987), ‘Turbulence statistics in fully de- 
veloped channel flow at low Reynolds number,” J. Fluid Mech., Vol. 177, pp. 
133-166. 

Keller, H.B. (1977), “Numerical Solution of Bifurcation and Nonlinear 
Eigenvalue Problems,” Applications of Bifurcation Theory, PH. Rabinowitz, 
ed., Academic Press, pp. 359-384. 

Lafon, A. and Yee, H.C. (1991), “Dynamical Approach Study of Spurious 
Steady-State Numerical Solutions for Nonlinear Differential Equations, Part 
III: The Effects of Nonlinear Source Terms in Reaction-Convection Equations,” 
NASA Technical Memorandum 103877, July 1991; International J. Comput. 
Fluid Dyn., Vol. 6, pp. 1-36, 1996. 

Lafon, A. and Yee, H.C. (1992), “Dynamical Approach Study of Spurious 
Steady-State Numerical Solutions of Nonlinear Differential Equations, Part IV: 
Stability vs. Numerical Treatment of Nonlinear Source Terms,” ONERA-CERT 


Building Blocks for Reliable Simulations 


37 


Technical Report DERAT 45/5005.38, also. International J. Comput. Fluid 
Dyn., Vol. 6, pp. 89-123, 1996. 

LeVeque, R.J. and Yee, H.C. (1990), “A Study of Numerical Methods for 
Hyperbolic Conservation Laws with Stiff Source Terms,” J. Comput. Phys., 
Vol. 86, pp. 187-210. 

Maday, Y. and Patera, A.T., (1989) “Spectral Element Methods for the 
Navier-Stokes Equations”, State of the Art Survey in Computational Mechanics, 
A.K. Noor, ed., ASME, New York, pp. 71-143. 

Maday, Y., Patera, T. and R0nquist, E.M., (1990) “An operator-integration- 
factor splitting method for time-dependent problems: Application to incom- 
pressible fluid flow, newblock J.Sci.Comput., Vol. 5, pp. 263-292. 

Minion, M.L. and Brown, D.L., (1997) “Performance of Under-resolved 
Two-Dimensional Incompressible Flow Simulations II”, J. Comput. Phys. 
138 , 734-765. 

Moin, P., and Kim, J. (1982), “Numerical investigation of turbulent channel 
flow,” J. Fluid Mech., Vol. 118, pp. 341-378. 

Moore, D.R., Weiss, N.O. and Wilkins, J.M. (1990), “The Reliability of 
Numerical Experiments: Transitions to Chaos in Thermosolutal Convection,” 
Nonlinearity, Vol. 3, pp. 997-1014. 

Moretti, G. and Abbett, M., (1966), “A Time-Dependent Computational 
Method for Blunt Body Flows,” AIAA Journal, Vol. 4, pp. 2136-2141. 

NEKTON User’s Guide, Version 2.8, 1991, Nektonics Inc., Cambridge, MA. 

Orszag, S. (1971), "Accurate solution of the Orr-Sommerfeld stability equa- 
tion," J. Fluid Mech., Vol. 50, pp. 689-703. 

Poliashenko, M. and Aidun, C.K. (1995), “Computational Dynamics of Or- 
dinary Differential Equations,” Intern. J. Bifurcation and Chaos, Vol. 5, pp. 
159-174. 

Schreiber, R. and Keller, H.B. (1983), “Spurious Solution in Driven Cavity 
Calculations,” J. Comput. Phys., Vol. 49, pp. 310-333. 

Shroff, G.M. and Keller, H.B. (1993), Stabilisation of unstable procedures: 
The RPM, SIAM J. of Numer. Analy., Vol. 30, No. 4, pp. 1099-1120. 

Shubin, G.R., Stephens, A.B. and Glaz, H.M. (1981), “Steady Shock Track- 
ing and Newton’s Method Applied to One-Dimensional Duct Flow,” J. Comput. 
Phys., Vol. 39, pp. 364-374. 

Sjogreen, B. and Yee, H.C., (2000) “Multiresolution Wavelet Based Adaptive 
Numerical Dissipation Control for Shock-Turbulence Computations”, RIACS 
Report 01.01, NASA Ames research center (Oct 2000). 

Sjogreen, B. and Yee, H.C., (2001) “On Entropy Splitting, Linear and Non- 
linear Numerical Dissipations and Long-Time Integrations, Proceedings of the 


38 


5th Intemat. Conf. on Spectral and High Order Methods, Uppsala, Sweden, 
June 1 1-15, 2001. 

Spalart, P.R., Moser, R.D. and Rogers, M.M. (1991), “Spectral methods for 
the Navier-Stokes equations with one infinite and two periodic directions,” J. 
Comput. Phys., Vol. 96, p297. 

Stephens, A.B. and Shubin, G.R. (1981) “Multiple Solutions and Bifurca- 
tion of Finite Difference Approximations to Some Steady Problems of Fluid 
Dynamics,” SIAM J. Sci. Statist Comput., Vol. 2, pp. 404-415. 

Sweby, PK. and Yee, H.C. (1994), “On the Dynamics of Some Grid Adapta- 
tion Schemes,” Proceedings of the 4th International Conference on Numerical 
Grid Generation in CFD and Related Fields, University College of Swansea, 
UK, also RIACS Technical Report 94.02, Feb. 1994. 

Sweby, P.K., Lafon, A. and Yee, H.C. (1995), “On the Dynamics of Com- 
puting a Chemically Relaxed Nonequilibrium Flow,” presented at the ICFD 
Conference on Numerical Methods for Fluid Dynamics, April 3-6, 1995, Ox- 
ford, UK. 

Thompson, J.M.T. and Stewart, H.B. (1986), Nonlinear Dynamics and 
Chaos , John Wiley, New York. 

Torczynski, J.R. (1993), “A Grid Refinement Study of Two-Dimensional 
Transient Flow Over a Backward-Facing Step Using a Spectral-Element 
Method,” FED- Vol. 149, Separated Flows, ASME 1993, J.C. Dutton and L.P 
Purtell, editors. 

Yee, H.C., Sweby, PK. and Griffiths, D.F. (1991), “Dynamical Approach 
Study of Spurious Steady-State Numerical Solutions for Nonlinear Differential 
Equations, Part I: The Dynamics of Time Discretizations and Its Implications 
for Algorithm Development in Computational Fluid Dynamics,” NASA TM- 
102820, April 1990; J. Comput. Phys., Vol. 97, pp. 249-310. 

Yee, H.C. and Sweby, P.K. (1992), “Dynamical Approach Study of Spurious 
Steady-State Numerical Solutions for Nonlinear Differential Equations, Part 
II: Global Asymptotic Behavior of Time Discretizations,” RNR-92-008, March 
1992, NASA Ames Research Center; also International J. Comput. Fluid 
Dyn., Vol. 4, pp. 219-283, 1995. 

Yee, H.C. and Sweby, P.K. (1993), “Global Asymptotic Behavior of Iterative 
Implicit Schemes,” RIACS Technical Report 93.11, December 1993, NASA 
Ames Research Center, also Intern. J. Bifurcation & Chaos, Vol. 4, pp. 
1579-1611. 

Yee, H.C. and Sweby, P.K. (1995), “On Super-Stable Implicit Methods and 
Time-Marching Approaches,” RIACS Technical Report 95.12, NASA Ames 
Research Center, July 1995; also, Proceedings of the Conference on Numerical 
Methods for Euler and Navier-Stokes Equations, Sept. 14-16, 1995, University 


Building Blocks for Reliable Simulations 


39 


of Montreal, Canada; International J. Comput. Fluid Dyn. Vol. 8, pp. 
265-286, 1997. 

Yee, H.C. and Sweby, P.K. (1996), “Some Aspects of Numerical Uncertain- 
ties in Time Marching to Steady-State Computations,” AIAA-96-2052, 27th 
AIAA Fluid Dynamics Conference, June 18-20, 1996, New Orleans, LA., 
AIAA J., Vol. 36, No. 5, pp. 712-724, 1998 

Yee, H.C., Torczynski, J.R., Morton, S.A., Visbal, M.R. and Sweby, RK. 
(1997), “On Spurious Behavior of CFD Simulations,” AIAA 97-1869, Pro- 
ceedings of the 13th AIAA Computational Fluid Dynamics Conference, June 
29 - July 2, 1997, Snowmass, CO.; also International J. Num. Meth. Fluids, 
Vol. 30, pp. 675-711, 1999. 

Yee, H.C. and Sweby, P.K. (1997), “Dynamics of Numerics & Spurious 
Behaviors in CFD Computations,” Keynote paper, 7th ISCFD Conference, 
Sept. 15-19, 1997, Beijing, China, RIACS Technical Report 97.06, June 1997. 

Yee, H.C., Sandham, N.D. and Djomehri, M.J., (1999) ‘ ‘Low Dissipative 
High Order Shock-Capturing Methods Using Characteristic-Based Filters”, J. 
Comput. Phys., 150 199-238. 

Yee, H.C., Vinokur, M. and Djomehri, M.J., (1999) “Entropy Splitting 
and Numerical Dissipation,” NASA Technical Memorandum 208793, August, 
1999, NASA Ames Research Center; J. Comput. Phys., 162, 33 (2000). 

Yee, H.C. and Sjogreen, B., (2001a) “Adaptive Numerical-Dissipation/Filter 
Controls for High Order Numerical Methods,” Proceedings of the 3rd Interna- 
tional Conference on DNS/LES, Arlington, Texas, August 4-9, 2001. 

Yee, H.C. and Sjogreen, B., (2001b) “Designing Adaptive Low Dissipative 
High Order Schemes for Long-Time Integrations,” Turbulent Flow Computa- 
tion, (Eds. D. Drikakis & B. Geurts), Kluwer Academic Publisher; also RIACS 
Technical Report, Dec. 2001 . 


